Multiscale Simulation of History Dependent Flow in Polymer Melt 
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We have developed a new multiscale simulation technique to investigate history-dependent flow 
behavior of entangled polymer melt, using a smoothed particle hydrodynamics simulation with 
microscopic simulators that account for the dynamics of entangled polymers acting on each fluid 
element. The multiscale simulation technique is applied to entangled polymer melt flow around 
a circular obstacle in a two-dimensional periodic system. It is found that the strain-rate history- 
dependent stress of the entangled polymer melt affects its flow behavior, and the memory in the 
stress causes nonlinear behavior even in the regions where Wi < 1. The spatial distribution of the 
entanglements {Z) is also investigated. The slightly low entanglement region is observed around the 
obstacle and is found to be broaden in the downstream region. 

PACS numbers: 83.10.Rs,83.10.Mj,83.60.Df,83.80.Sg 
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Industrial products using polymeric materials have be- 
come increasingly integral to our lives. One of the impor- 
tant characteristics of thermo-plastic polymeric material 
is that it can be easily molded and processed by control- 
ling its state from a solid to a melt, which is beneficial 
in a variety of practical applications. The melt state 
of polymeric material can exhibit a lot of characteristic 
phenomena, e.g., die-swell and rod-climbing[lj, depend- 
ing on the dynamic response of the polymer's microscopic 
internal states under an imposed strain or strain-rate. 

Predicting the flow of polymeric fluid is difficult be- 
cause the macroscopic flow behavior depends heavily on 
the dynamic behavior of the microscopic internal states, 
and the states of polymers with very high molecular 
weight are complex at the microscopic level. When 
a polymer melt consists of polymers with a molecular 
weight M higher than the entanglement molecular weight 
Me, the polymer molecules are entangled with each other, 
and the relaxation time of the polymer conformation is 
long compared to that of dilute polymer solutions be- 
cause of the entanglement. Therefore, it is difficult to 
predict the rheological behavior, even for a homogeneous 
bulk entangled polymer melt, using molecular dynamics 
approaches without any coarse-graining procedures for 
the microscopic internal degrees of freedom. Reptation 
theory [i-!3|, in which a polymer chain is coarse-grained to 
a tube confined by the surrounding polymer chains, ex- 
plains the complex behavior of entangled polymer melt 
composed of mono-disperse linear polymers. To apply 
the concept of a tube to a wider variety of polymers, 
the original reptation theory has been improved by in- 
troducing several important physical mechanisms [5|-l8|. 
Extended reptation theory based on the Fokker-Planck 
equation for tube segments succeeds in explaining many 
experimental resultsfS-'s'] . However, it is difficult to ap- 
ply the theory to arbitrary polymer melts because it is 
too mathematically complex to incorporate the molecular 
pictures of polymers of arbitrary architectures into the 



Fokker-Planck equation. On the other hand, Langevin 
stochastic models based on reptation theory using the 
complex molecular configurations have been developed 
and can reproduce the rheological properties of vari- 
ous polymeric materials with branches and/or molecular 
weight distributions 'q'-T^] . These models use efficient nu- 
merical computation to predict the bulk rheology of poly- 
mer melts; however, they are not applicable to macro- 
scopically inhomogeneous flows only by themselves. The 
macroscopic flow behavior of polymer melt is usually pre- 
dicted using a fluid dynamics simulation with a consti- 
tutive equation that describes a nonlinear and history 
dependent relationship between the stress and strain- 
rate flelds to represent complex fluids without micro- 
scopic details. However, no general constitutive equation 
that is applicable to entangled polymer melts with arbi- 
trary polymer architecture exists because a general poly- 
mer melt, e.g., polymers with branches and/or molecu- 
lar weight distributions, has extremely complex molecu- 
lar states, which makes the stress response unpredictable 
even in simple flow patterns. As described above, each of 
the microscopic and macroscopic approaches has limited 
applications. Recently, a molecularly derived constitu- 
tive equation technique is developed for low-molecular 
polymer melts without entanglements [13]. If it is gener- 
alized to a variety of polymer architectures, the technique 
can be attractive and useful. However, it is in an early 
phase of development toward generalization. 

As an alternative to using a constitutive equation, we 
propose a new simulation technique that incorporates 
Langevin course-grained simulators [9l-fl2| into a fluid dy- 
namics simulation. This approach is a type of micro- 
macro simulations that were mainly developed for dilute 
polymer solutions in which the stress tensor of each fluid 
element is obtained from a microscopic simulation [ij- 
and applied not to dilute polymer solutions but to 
well entangled polymer melts here. In contrast to the 
dilute polymer solutions, the polymer melts have en- 
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tanglement interactions between neighboring polymers. 
In order to investigate the entangled polymer melts us- 
ing the micro-macro simulation, we need to satisfy (i) 
a local stress tensor is represented by an ensemble of a 
sufhciently large number of entangled polymers and to 
consider (ii) strong correlation among the polymers in 
an ensemble because of the entanglements. The advec- 
tion of the stress, (v • V)cr, in the Eulerian description 
is well established in any fields. However, it only ad- 
vects the macroscopic quantity of the stress (the statis- 
tical conformation tensor) without involving microscopic 
molecules themselves. Because the correlations among 
polymers, namely the entanglements, intrinsically affect 
their macroscopic properties in the entangled polymer 
melts, the advection of an ensemble of polymers cor- 
relating at the microscopic level should be consistently 
managed with the advection of the macroscopic vari- 
ables. Therefore, we adopt a Lagrangian fluid particle 
method TMi to ensure the advection of the ensemble 
of entangled polymers. Each ensemble of entangled poly- 
mers is assigned on each fluid element. The advection of 
microscopic details is essential to the description of the 
macroscopic flow behavior for entangled polymer melts 
because the conformations of entangled polymers are con- 
strained by and correlated to the surrounding polymers. 
Note that the Eulerian techniques can still be useful even 
for entangled polymer melts when the system has a trans- 
lational symmetry [2C' 
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To maintain the information pertaining to entangle- 
ment and deformation in polymers, we perform a La- 
grangian fluid particle simulation in which each fluid par- 
ticle has a microscopic level simulator that accounts for 
the internal states of the fluid particle 19j. Assuming 
that the polymer melt is an isothermal and incompress- 
ible fluid, the governing equations for the i-th fluid par- 
ticle that constitutes the polymer melt are given by the 
following equations; 

P.^ = v•<T,-Vp, + ^^ ^ = v„ (1) 



dt 



dt 



(2) 



where F'' is a body force and the pressure Pi = Pi{{pi}) 
is properly considered to guarantee the incompressibility. 
These equations are solved using macroscopic variables, 
except for Eq. ^ . The local stress tensor CTi depends on 
a microscopic ensemble Qi of entangled polymers, which 
represent entangled states of polymers under an instan- 
taneous local strain-rate Hi = (Vv)^. 

The slip- link model [HI is a simulation model that 
can accurately describe the dynamics of entangled poly- 
mers. The model is composed of conflning tubes with 
some entanglement points, called slip-links, which con- 
flne a pair of polymers and represent effective constraints 
in virtual space. The average number of slip-links or en- 
tanglements on a polymer at the equilibrium state is rep- 
resented as {Z)cq = M/Mc- In the simulation, we trace 



the conflgurations of conflning tubes constrained by the 
slip-links. The slip-links are relatively convected each 
other and the conflning tubes are deformed according to 
the macroscopically obtained local velocity gradient ten- 
sor n. The reptations of polymers generate or eliminate 
slip-links. For given chain configurations, the stress ten- 
sor crP derived from deformations of entangled polymers 
is calculated from cr^^ = ^^('"ja'^^/s/k^D/as where Os 
is the unit length of the slip-link model and r^^^ is the 
a-component of the j-th tube segment vector connecting 
adjacent slip-links along a polymer. The unit of stress do 
in the slip-link model relates to the plateau modulus Gn 
as follows: = {15/4)GN[m- The slip-link model has 
two characteristic time-scales: the Rouse relaxation time 
and the longest relaxation time tj. The Rouse relax- 
ation time tr and the longest relaxation time relate to 
(Z)eq as follows: tr = {Z)l^U and Td oc (Z)3q'^^[|^ [JJ], 
where to is the time unit of the slip-link model. The con- 
tour length relaxation of a conflning tube occurs on the 
time-scale of tt{ , while the orientational relaxation occurs 
on the time-scale of Td. These two characteristic times 
appear in the stress relaxation. 

Each polymer simulator describing a fluid particle com- 
putes the polymer conflgurations at each time step, and 
the recorded conflgurations are used as the initial con- 
ditions of the next time step. Typically, the macro- 
scopic time unit tmacro and microscopic time unit <micro 
have a large time-scale gap, and therefore the macro- 
scopic time unit tmacro must be divided into iVtmicro- Be- 
cause the slip-link model used here is sufficiently coarse- 
grained and the time unit te can be the same as the 
time-scale of the macroscopic fluid tmacro, we employ 
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Note that we set the local stress of the macroscopic 
fluid to <T = o-P -I- <t'^, where cr'^ is an extra dissipative 
stress tensor. Because the slip-link model is a Langevin 
coarse-grained model based on reptation theory, micro- 
scopic dynamics smaller than a tube segment are treated 
as a random force exerted on a slip-link, and the contri- 
bution from the microscopic dynamics does not explic- 
itly appear in the stress tensor of the slip-link model. 
We assume the dissipative stress cr'^ to be the Newtonian 
viscosity ?7jD, where D = /t + k?- is a strain-rate tensor. 

The main procedures of our simulation are summa- 
rized as follows: (1) Update {v^}, {vi] at the macro- 
scopic level. (2) Calculate {pi]^{Ki] at the macroscopic 
level. (3) Obtain {pi\ from the density distribution {pi\ 
at the macroscopic level. (4) Update the local slip-links 
of the z-th fluid particle under the local strain-rate Hi 
and then obtain cTj from the resulting configuration of 
the slip-link model. This procedure is executed on each 
fluid particle in turn. (5) Calculate {V • ctj} and {Vp^} 
at the macroscopic level. (6) Return to (1). 

We update {v^}, {r^} by integrating Eq. H]). We 
calculate the density at the position of each particle 
in the new conflguration using a method in the usual 
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smoothed particle hydrodynamics technique [22|: pi = 
^j'nijW{\Yj — Yi\,h), where is the mass of the i- 
th particle and W^(|r|,/i) is a Gaussian-shaped function 
with cutoff length 2h. The deviation of the local den- 
sity from the initial constant density po results in a lo- 
cal pressure force — Vp. To obtain the spatial deriva- 
tive of the velocity field, stress field, and pressure field 
(Vv, V • cr, Vp), we use a technique that was developed 
for modified smoothed particle hydrodynamics [23l |2J|. 

To demonstrate the efficiency of the proposed multi- 
scale simulation, we consider a system in which the flow 
history can affect the flow behavior. One such system 
is a polymer melt flow around an infinitely long cylinder 
oriented in the z-direction with a radius Tc, which flows 
in the x-direction. Because of the symmetry of the sys- 
tem, we can treat the system as two-dimensional, and 
the flow can be described as two-dimensional flow in the 
xy-plane. We assume a non-slip boundary condition for 
the velocity on the surface of the cylinder and periodic 
boundary conditions at the boundaries of the system. 
The dimensionless parameters governing the problem are 
the Reynolds number Re = p\v\rc/rf'^ the Weissenberg 
number Wi = TdD^y, and the viscosity ratio rjd/rj^, where 
V is the average flow velocity. The zero shear viscosity 
r]'^ of the polymer melt is given by ?7p + ri^, where 7]^ is 
the zero shear viscosity of a polymer melt described only 
by the slip-link model. From the rheological data shown 
in Fig. [T] (a) , which were derived from the slip- link simu- 
lation with {Z)eq = 7 in the bulk, we obtain the longest 
relaxation time ~ 200te and the zero shear viscos- 
ity ~ 17.5aete- The cylinder radius Tc was set to 3a, 
where a is the unit length in the fluid particle simulation, 
and we assign the unit mass m to all fluid particles. The 
wall of the cylinder consists of fixed fluid particles evenly 
spaced on the perimeter. Each fluid particle consists of 
f 0000 polymers, enough to describe the bulk rheological 
properties of the polymer melt under an imposed shear 
and/or extensional deformation [1^]. About 900 fluid par- 
ticles with 10000 polymers each are evenly placed initially 
in the system and then move according to Eq. ([!]); there- 
fore, we need to simultaneously solve the dynamics of 
9000000 polymers. Because the diffusive motion of the 
center of mass of a single polymer is negligible compared 
to the translational motion of the center of mass of an 
ensemble of entangled polymers, we neglect any trans- 
portation of polymers between adjacent fluid particles. 
With this assumption, each slip-link simulation can be 
performed independently of the others, making parallel 
computing effective in this multiscale simulation. 

Under the body force / {tj'^ / aQ = (5.0 x 10-^0), 
the flow becomes steady-state in about lOOOtc- The aver- 
age flow velocity v in steady-state in this system is nearly 
equal to (0.04, 0)a/te for a fully Newtonian flow 77° = ij^ 
(o-P = 0), and (0.055, 0)a/tc for a polymer melt flow with 
{Z)eq = 7 and 'qd/ff = 0.1. In both cases. Re is less than 
0.2, and the flow is laminar. In the polymer melt case. 




FIG. 1. (a) Rheological data and (b) relaxation modulus G{t) 
obtained from the slip-link simulation with (Z)eq = 7 in the 
bulk. The storage modulus G' (dashed line), the loss mod- 
ulus G" (dotted line), and the complex viscosity r)* (solid 
line) are plotted against the angular frequency uj in (a). The 
magnitude of the complex viscosity, |J7*| = yjG'^ 4- C'^/oj, is 
considered to be the shear viscosity r] of the low shear-rate 
region (Cox-Merz rule). The normalized relaxation modulus 
G{t)/G{0) (solid line) and a single exponential curve with a 
relaxation time r (dashed line) are plotted against the nor- 
malized time t/r in (b). 
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FIG. 2. Color contour maps of time-averaged magnitudes of 
the velocity field |v|[a/te], strain-rate [1/te], and shear 
stress over the zero shear viscosity liyxyl/ri^i^/te] for (a) the 
polymer melt with {Z)cq ~ 7 and Jyd/j?" = 0.1 and (b) the 
Newtonian fluid with r;" = r^d in steady-state. Figure (c) 
shows (Z) in the polymer melt. The regions inside the red 
lines in (a) correspond to Wi > 1 and those in (c) to rnj > 1 
where 7 = ^^/Tt^k^ k). 



the average flow velocity is higher than that of Newto- 
nian flow, i.e., the flow exhibits shear thinning behavior 
because of Wi > 1 in the vicinity of the cylinder. 

To investigate the velocity field v, strain-rate D^y, and 
shear stress axy in steady-state, we employ a linear inter- 
polation to transform the data at the particle positions 
into the values at regular lattice points and then time- 
average the data evaluated at the lattice points. To elimi- 
nate the noise of the data, the time-averaging was carried 
out from lOOOic to 2000io- Figure [2] shows the spatial dis- 
tributions of |v|, \Da:y\ and IcTxyl/rj^ in steady-state for 
(a) the polymer melt with {Z)cq = 7 and rj^/rf = 0.1 
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FIG. 3. Stream lines (solid line(A) and dashed line(B)) ob- 
tained from the time-averaged velocity field with particle 
tracking are plotted against position in (a). These stream 
lines are in the region where Wi < 1. The strain-rate D^y 
and the shear stress a:cy along the stream lines are plotted 
against the elapsed particle tracking time in (b). Each line 
type corresponds to each stream line. 



and (b) the Newtonian fluid with 77° = ?yd- 

Reflecting laminar behavior, the magnitudes of v and 
in Fig. [2] appear to be nearly symmetric between 
the upstream and downstream regions unlike In 
Fig- 12] (a): the nonlinear relationship between Oxy and 
Dxy is observed near the cylinder where Wi > 1 because 
of shear thinning. Moreover, the shear stress a^y of the 
polymer melt exhibits an apparent asymmetry between 
the upstream and downstream regions. In general, the 
viscoelastic stress relates to the shear-rate through the 
relaxation modulus according to the following equation: 
a(t) = J_^dt'G{t — t')j{t'); the stress depends on the 
flow history. We investigate the behavior of (Txy and D^y 
on two typical stream lines (A: solid line and B: broken 
line shown in Fig. [3] (a)) in the region where Wi < 1. 
Figure [3] (b) shows axy and Dxy along these stream lines 
and plots them against the elapsed particle tracking time. 
The time origin t = is set to tmin, when the strain-rate 
Dxy is at a minimum in each stream line. A nonlin- 
ear relationship between axy and Dxy is not evident in 
either stream line because Wi < 1; however, the mini- 
mum of axy appears to be shifted from that of Dxy by a 
time difference Ia- Applying a step shear strain 7 = 0.5 
to the slip-link model with (Z)eq = 7 in the bulk, we 
obtain the relaxation modulus G{t) = a{t)/^ shown in 
Fig. [T] (b). The stress relaxation time r is estimated to 
be r ~ 50ic(=i tr) when G{t = t)/G(0) = 1/e. The time 
difference is found to be nearly equal to r. 

Finally, we investigate the spatial distribution of en- 
tanglements (Z) shown in Fig. [2] (c). In the vicinity of 
the cylinder, the entanglements slightly decrease. When 
trJ > 1, polymer's contour length is highly extended. 
The extended polymer longer than {Z)cqa is easy to 
shrink, which causes disentanglement. In the upstream 
and downstream regions around the cylinder, the reduc- 
tion of (Z) is also observed because of nonzero \Dyy\. The 
flow advection broadens the region where {Z) < {Z)cq in 
the downstream region. The tail length in the down- 
stream region is roughly estimated to be tr|v| ~ rc. 



Using the Lagrangian particle method to trace and 
maintain the entire configurations of polymers, we have 
been able to describe the memory effect in the polymer 
melt flow. The presented multiscale simulation is ap- 
plicable to various polymer melts, because the slip- link 
model or the alternative course-grained models can 
address a variety of polymer architectures, e.g., linear 
and/or branched polymers, polymer blends, and polydis- 
persed polymers. The multiscale simulation is advanta- 
geous because it employs a fully Lagrangian method at 
the macroscopic level, while conventional micro-macro 
techniques which have difficulties accounting for the 
macroscopic advection of microscopic internal states. 
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